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The polarizability measures how the system responds to an applied electrical field. Computa- 
tionally, there are many different ways to evaluate this tensorial quantity, some of which rely on 
the explicit use of the external perturbation and require several individual calculations to obtain 
the full tensor. In this work, we present some considerations about symmetry that allow us to take 
full advantage of Neumann's principle and decrease the number of calculations required by these 
methods. We illustrate the approach with two examples, the use of the symmetries in real space 
and in spin space in the calculation of the electrical or the spin response. 



I. INTRODUCTION 

The redistribution of electrons in a finite system that 
occurs when it is exposed to an external electromag- 
netic field is characterized by a set of constants called 
polarizabilitiesi. The static polarizabilitics refer to static 
fields, whereas the dynamical polarizabilitics arc func- 
tions of the frequency of the external field. Usually the 
name polarizability is restricted to the constants that de- 
termine the cited redistribution to first order in the ap- 
plied field - when high-intensity fields are applied, one 
needs to make use of the higher-order ones or hyperpo- 
larizabilities. 

Besides characterizing the electromagnetic field-matter 
interaction, they arc also important when studying colli- 
sion phenomena, and as coarse indicators of physical size, 
structure and shape. The knowledge of the polarizability 
helps to obtain numerous physical quantities that depend 
on it: the dielectric constant and the refractive index of 
extended systems, the long-range interaction potentials 
between polarizable systems, van der Waals constants, 
etc. The dynamical polarizability, also, contains precious 
information about the excited states of a finite system: 
the peaks of the polarizability as a function of the energy 
determine the excitation energies of the system. The ab- 
sorption cross section of a finite system is trivially related 
to the imaginary part of the dynamical polarizabilityii^. 
The oscillator strengths of the transition are the weights 
of the peaks of this absorption spectrum. 

The polarizabilities may also be classified as dipole, 
quadrupole, etc., depending on the shape of the perturb- 
ing field: For example, the dipole polarizabilitics charac- 
terize the response of the system to a dipole field. Rig- 
orously speaking, one should also specify the physical 
magnitude that is measured, and speak, for example, of 
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the dipolc-dipolc polarizability: the response of the sys- 
tem dipole to a dipole field. In this Article, we will ex- 
clusively deal with dipolc-dipole first-order polarizabil- 
ities - the most commonly studied ones, and arguably 
the most important for nanostructure characterisation. 
Since there are three possible directions for the perturb- 
ing dipole field, and three possible directions for the sys- 
tem dipole response, the dipolc-dipolc polarizability is in 
fact a three-dimensional tensor. 

Because of its tensorial character, one way of simpli- 
fying the calculation of the polarizability is to use Neu- 
mann's principle^.: the polarizability tensor of the system 
must be left invariant under any transformation that is 
also a point symmetry operation of the system. This 
condition of invariance reduces the number of indepen- 
dent tensor components, since it signifies relationships 
between those components, thus potentially reducing the 
number of calculations necessary to obtain the full tensor. 

Numerous theoretical techniques can be used to calcu- 
late polarizabilitics, with varying level of accuracy and 
detail. In particular, there is a class of methods that rely 
on the explicit use of the external perturbation, i.e. each 
line of the tensor is obtained by performing one calcu- 
lation. Because of this, when using these methods it is 
not always obvious along which directions the perturbing 
fields should be aplied in order to make full use of Neu- 
mann's principle. In this Article we discuss how this can 
be done. 

Stricto sensu, the polarizabilities refer to the reac- 
tion to spin- independent (i.e., electrical) perturbations 
measured by spin-independent obscrvablcs; they are re- 
ferred to as density-density response functions. However, 
one can also think of more general response functions 
and apply spin-dependent perturbations and/or look at 
spin-dependent obscrvablcs, obtaining in this way spin- 
density, density-spin and spin-spin response functions - 
these generalized objects are sometimes referred to as 
susceptibilities. Even though this work is mainly con- 
cerned with the polarizability itself, we will consider these 
more general objects, since we will also show how the cal- 



culation of the density-density and spin-spin response of 
a spin-saturated molecular system may be simplified. 

In the following section, we recall the necessary defini- 
tions, and list some of the first-principles techniques that 
can be used to calculate polarizabilities of tcchnologicaly 
relevant nano and bio structures - some of which can 
benefit from the simplifications proposed in this Article. 
These are discussed in Section Hill Section HVl shows how 
essentially the sames ideas may serve to simplify the cal- 
culation of the singlet and triplet excitations of systems 
whose ground-state is spin-saturated. Finally, we present 
two systems where the technique was used to compute 
the density-density and the spin-spin responses. 



II. THE POLARIZABILITY TENSOR 

We will now introduce the notation and key quantities 
that are relevant for the use of symmetries to obtain the 
difi'erent linear response functions. Let Ug- (cr ^tii) be 
the spin-density of a system of N electrons: 

N 

= mt)\J2S{r-h)6ggjm), (1) 

If we apply an infinitesimal perturbation, 5vu(f, oj) (a =] 
,J,), the response of the density, 5na(f,t) = na{r^t) — 
ncr(r, 0), will be related to the perturbation, to first or- 
der, by a density response function, x- In tlie frequency 
domain this is expressed as; 

5n,{r,u:) = Y, /dV x..' (r, r', c.)<5«,, (r', c.) . (2) 



The variation of the total density^ n ~ n| + n^, and of 
the magnetization density, m = n-\ — n^, are given by: 

(5n(f, (jj) = (5n|(r, a;) H- (5n|(r, cj) , (3a) 
(5m(f, cj) = (5n|(r, a;) — (5n|(r, cj) . (3b) 

After the system is perturbed, one can obtain infor- 
mation about the system by looking at the variation of 
some observable: 5{d){t) = (d)}(t) - (d))(0). In our case, 
we will be looking at the dipolc of the system in each of 
the spatial directions, Xi. In order to learn about the 
spin modes, one may also look at A^ct^, where ct^ is the 
Pauli z-matrix. In the frequency domain, the behavior 
of these observables is given by: 



5{Xi){uj)= J d^rXiSn{r,uj) . 
S{Xi(Tz){uj) — d^r XiSm{r,uj) 



(4a) 
(4b) 



One also has to define which kind of perturbations are 
to be considered: We will restrict hereafter the formula- 
tion to dipole perturbations of two kinds: 
(i) Spin-independent perturbations of the form: 



In this case, the variations Sn and Sm are: 

Jn(f,w) = -k(w) Jd^r'x^"''\f,f',u;)x'j, (6a) 

Sm{r,uj) = -K(tj) y"dV x[™"l(r,7=',w)a;^. , (6b) 

where we have defined the new objects: 

x'""' =XTT+Xn+XiT+XU (7a) 

x'™"' =XTT+XTi-XiT-XU (^b) 



The observables 6{Xi){uj) and S{Xi(Tz){uj) will be given 



by: 

S{X,)^;\lj) = -niiv) J l'd^rd\'x,x^^"^ry,iu)x'^ , 

(8a) 

5{X,az)f\uj) - -^Kico) J jd\d\' x,x^'"'^\r,7^,uj)x'^ . 

(8b) 

The superscript "[n]" means that the observable is mea- 
sured after a spin-independent perturbation of the form 
given in Eq. ([5]) has been applied, whereas the subscript 
"j" means that this perturbation has been applied in the 
direction j. 

(ii) Spin- dependent perturbations of the form: 

Jz;[™l(^^) = | "''■'■1''! ' (9) 
[ XjK(u;) , (J =1 . 

Or, making use of the Pauli z-matrix: 

6v^"'\f,uj)^-XjK{Lj)az. (10) 

The variations of n{f,uj) and m{f,uj) are: 

Sn{r,u;) = ^k{lj) ^dV x'""' (r^ , w)^ , (11a) 
Sm{f,iu) = y"dVx[™™l(^^,w)x; , (lib) 

with the definitions: 

x'""' = XTT - xn + xn - XU (12a) 



X' 



= XTT - XTi. - XJ.T + Xii (12b) 



Now, the observables S{Xi){u!) and 6{Xi(Jz){uj) will be 
given by 



S{X,piu;) = -Kiu;) jj. 



d^rd^rx,x^''"'\r.r',uj)x' 



5{X,azYp{u:) = ^k{lo) j j d\d\x,x^"'-'\r,7^ ,lo)x'^ 



5v^:\r,Lu)^^XjK{u). 



(5) 



(13a) 
(13b) 
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We now look at the quotients between the induced vari- 
ations 5{Xi){Lu) and S{Xiaz){(jj) and the strength of the 
perturbation, k{uj), for each one of the cases: 



a) 



k{ll!) 

Sjx.a^pju;) 



(14a) 
(14b) 
(14c) 
(14d) 



The usual definition of the polarizability refers to the first 
of these expressions, aij = a,[""' . However, one may also 
be interested in the other kinds of responses. We will 
use the same notation (a) for these response functions, 
although the name polarizability should be restricted for 
the first case. Another way of defining these functions is: 



(15) 



Obviously, the two definitions are related, and one may 



retrieve the a: 



components from the and vice- 
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Regarding its spatial structure - and thus dropping 
the spin indexes - the dynamical polarizability elements 
may be arranged in a second-rank 3x3 symmetric tensor 
a{LL!). The absorption cross-section tensor is then simply 
proportional to its imaginary part: 



cr{uj) — Jma(uj) , 



(17) 



where c is the speed of light. 

The static or the dynamical (hyper)polarizabilities can 
be calculated within different theoretical approaches. 
Our main concern here is (time-dependent) density func- 
tional theory - (TD)DFT— , but our arguments are 
quite general and apply equally well to other electronic 
structure methods. The simplest of these is perhaps 
obtaining the static (hyper)polarizabilities through fi- 
nite differences, i.e., by applying small static electrical 
fields Ej, and then calculating numerically the deriva- 
tives. On the other hand, the dynamical polarizabil- 
ity can be calculated through real-time propagation of 
the time-dependent Kohn-Sham cquations^i. Moreover, 
both static and dynamical polarizabilities can be ob- 
tained through straighforward perturbation theory. In 



this case, "perturbation theory" refers to the applica- 
tion of Sternheimer equation^i^ in one way or another, 
either for the static or for the dynamical^i^ case. Another 
very recent and quite promising approach, is a efficient 
Lanczos-based method^. 

For all these methods the calculations rely on the ex- 
plicit use of the external perturbation. If we want to 
calculate the full tensor, we have to perform three calcu- 
lations, one for each spatial direction. If, moreover, we 
need both the density and spin modes, we have to make 
two calculations per spatial direction, one for each of the 
two possible perturbations discussed above. In the next 
sections we will show how the number of actual calcula- 
tions can be severely reduced when using these methods 
by taking advantage of the symmetries of the system. 

Note that the polarizability can also be obtained, for 
example, from the response functions, using, for example, 
the formalism developed by M. Casida— liii^^iii^. This is a 
case in which the considerations presented below will not 
simplify the calculations, since they do not proceed by 
the successive application of the different perturbations. 



III. SPATIAL SYMMETRIES 

In this section, we will drop the spin indices, since the 
whole argument is valid for any of the spin components 
of the polarizability. 

By making use of Eq. ([15]), it is easy to prove that 
a{uj) is a proper tensor: If we consider a second or- 
thonormal reference frame {e'l, 62, e'^}, a{Lu) transforms 
following the tensorial transformation law: 



a'{uj) = P*a(w)P, 



(18) 



where Q;'(a;) is the polarizability in the second reference 
frame, and P is the rotation matrix between the two 
frames. 

This tensorial character of the polarizability permits 
us to work in any orthonormal reference frame; once we 
obtain its values, we may easily transform it by straight- 
forward matrix manipulation. We may then choose the 
frame which is most appropriate, bearing in mind the 
geometry of the molecule, and this can reduce the total 
number of calculations. However, this liberty does not 
allow us to make full use of symmetry. For this purpose, 
we need to work with non-orthonormal directions. 

Let us consider three linearly- independent, but possi- 
bly non-orthogonal, unit vectors {pi,P2,P3}- We define 
the polarizability elements Q!„t,(w) as: 



d-'rd'^r' {f-Pu)xiry,u;)ir' ■ p^) . (19) 



This corresponds to a process in which the polarization 
of the perturbing field is along pu, and the dipole is mea- 
sured along pu- If we know the 3x3 matrix a.{uj)^ we get 
the real tensor ol{u}) by making use of the following sim- 
ple relationship, which can be obtained once again from 
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Eq. (Uni): 



a{Lu) = P*q;(w)P . 



(20) 



P is the transformation matrix between the original or- 
thonormal reference frame and {pi,p2,P3}. Note that 
this transformation is in general not a rotation; P is 
not unitary. Moreover, no matter how familiar it looks, 
Eq. ((20|) is not a change of coordinates: (i(ijj) is not the 
polarizability tensor in the new reference frame. And 
finally, also note that the traces of a and a do not coin- 
cide: 



Tr [a{uj)] ^ Tr [P*a(w)P] = Tr [a(tj)PP*] 



(21) 



but PP* 1. Notwithstanding this, it tells us that we 
may obtain the polarizability tensor by calculating the 
related object d(w). 

Now let us assume that the molecule under study pos- 
sesses some non-trivial symmetry transformations - to 
start with, wc consider that it has two, A and B. We 
consider an initial unit vector, pi, and define: 



P2 = Api 

P3 = Bp2 



(22) 



We assume that this may be done such that the set 
{pi,P2,P3} is linearly independent. 

We then perform a TDDFT calculation with the per- 
turbing field polarized in the direction pi . This permits 
us to obtain the row {an, ai2-, 0:13}. Since the matrix is 
symmetric, we also have the column {fin, a2i, aai}. The 
symmetry of the molecule also permits us to obtain the 
diagonal: {533 = 022 = ctii}- The only missing element 
is 523 = C(32, but it is easy to prove that: 



"23 = Q;i,^-ir 



(23) 



which we can also get from our original calculation. The 
conclusion is that we have access to the full tensor by 
making only one calculation. 

To fix ideas, we use the example of a molecule with 
one n-th order axis of symmetry (n > 2). Let TZ be 
the rotation of 27r/ri degrees around this axis. We then 
choose pi not coUinear with this axis, and also not per- 
pendicular to it. If we define TZpi = p2 and TZp2 = Ps, 
the set {pi,P2,'P3} will be linearly independent. In this 
case, moreover, since A = B = TZ, Eq. (^5)) reduces to 
Q!23 = ai2. 

It may very well be that wc may only find two linearly 
independent "equivalent axis" , pi and p2 , related by a 
symmetry transformation, A - this is the case of a system 
that possesses only a plane of symmetry, or only an axis 
of symmetry of order two. We may then define pi to be 
a tilted vector with respect to this plane (not contained 
in it, and not perpendicular). Then, p2 = Api, where 
A is the reflection on the plane, is an equivalent vector, 
and P3 can be chosen to lie in the symmetry plane (the 
obvious choice will be P3 = pi /\ p2, that ensures linear 
independence). We then only need two calculations, one 



with the polarization along pi (or P2) and another with 
the polarization along ^3. Moreover, if pi is chosen to be 
tilted exactly ■n/2 with respect to the symmetry plane, 
the system of vectors is orthonormal, and we do not even 
need to apply Eq. ([20]) . Note that this case applies to all 
planar molecules. 



IV. SINGLET AND TRIPLET EXCITATIONS IN 
SPIN-SATURATED SYSTEMS 

We consider a system whose ground state is spin- 
saturated. It verifies: 



TT 
y, ■ 



And; in consequence, 



0, and 



[mm] _ 2q,TT 



2a 
2a 



n 



n 



(24a) 
(24b) 



(25a) 
(25b) 



Despite this symmetry, in order to obtain all spin- 
components (which in this spin-saturated case are only 
two independent ones), by making use of the two types 
of perturbations defined in Eqs. ([5]) and ([9]), we would 
still need two calculations: one perturbing with a spin- 
independent potential - in order to obtain , and one 

with a spin-dependent one - in order to obtain aj™"*^ . 

However, one can easily use a similar scheme to the 
one outlined in the previous section in order to calculate 
the two polarizabilities in only one shot. The idea is to 
apply a perturbation in the form: 



or, in the Pauli matrix language: 



(26) 



(27) 



It is then easy to verify that the response of the dipole 
observables will then be given by: 



5(X,)m(c.) = -i.(^)al7l 



S{X.a.)^P{u;) = -U-)a^"'^ 



(28a) 
(28b) 



thus providing us with the components of the two re- 
sponse functions that we need with only one calculation. 



V. EXAMPLES 

There are many complex molecules of tecnological rele- 
vance that present symmetries such that the schemes out- 
lined in the previous sections can be used. We chose two 




FIG. 1: (Color online) Two different views of protonated 
triphenylguanidine (point group Csh)- In this picture, the 
dark blue atoms are nitrogen, the cyan represent carbon, and 
finally the white atoms are hydrogen. 



of them to illustrated the method: protonated triphenyl- 
guanidine and one hidrogenated silicon cluster SigHig. 
Triphenylguanidine compounds are regarded as interest- 
ing for quadratic nonlinear optical applications while 
hidrogenated silicon is an important optico-electronic 
material with potentially important technological appli- 
cations. 

The ground-state of protonated triphenylguanidine is 
spin-saturated, has one proper axis of symmetry of order 
three, one plane of symmetry, and one improper axis of 
rotation of order three (see Fig. [T]). The ground-state 
of SisHis is spin-saturated, has an inversion center, one 
proper axis of symmetry of order three, three proper axis 
of symmetry of order two, three planes of symmetry, and 
one improper axis of symmetry of order six (see Fig. [2]). 
This means one can make use of the schemes outlined in 
the previous sections to obtain all the components of the 
response functions with only one calculation: all compo- 
nents of both a[""l(w) and a[™™l(cij) tensors. 

Without using the symmetry of the system the re- 
sponse functions were computed by applying spin- 
independent and spin-dependent perturbations from 
Eqs. ([5]) and ^ with polarization directions along the 
X, y and z directions. This way the response functions 
were straighforwardly obtained but required a total of six 



FIG. 2: (Color online) One view of SigHis (point group Dzd)- 
In this picture, the dark brown atoms are silicon and the grey 
atoms are hydrogen. 



time-propagations. 

To use the symmetry a set {pi,p2,P3} is needed. In 
the case of protonated triphenylguanidine we built it by 
defining pi to be a vector tilted 7r/4 with respect to the 
plane of symmetry and the two symmetry transforma- 
tions A and B to be an inversion with respect to the 
plane and a 27r/3 rotation around the axis of symmetry 
of order 3. In the case of SisHig we chose pi to be a 
vector tilted 7r/4 with respect to the axis of symmetry of 
order three and both symmetry transformations A and 
B to be 27r/3 rotations around the same axis. 

Applying a perturbation of the same form as Eq. (j26p 
with a polarization direction along pi and using Eqs. (j20p . 
pSap and (j28bp allowed us to obtain the response func- 
tions with just one calculation in both cases. 

All response calculations were done with the code 
octopu9i^ using the Perdew-Zunger— parametrization 
of the adiabatic local density approximation for the 
exchange-correlation potential. This method has al- 
ready been successfully used for the calculation of op- 
tical spectra in a variety of systems, ranging from small 
moleculesi^ and clustersii to biological systemsi^. 

To represent the wave- functions in real space we used 
a uniform grid with a spacing of 0.195 A and a box com- 
posed by spheres of radius 5 A around every atom. In 
order to propagate the Kohn-Sham orbitals we employed 
state of the art algorithmsiS,. A time step of 0.0048 fs 
assured the stability of the time propagation, and a total 
propagation time of 19.35 fs allowed a resolution of about 
0.1 eV in the resulting spectrum. 

The results obtained are summarized in Fig. [3] and 
Fig. m where we plot the averaged absorption cross- 
section (trace of the tensor) for singlets and triplets. 
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FIG. 3: (Color online) Protonated triphenylguanidine aver- 
aged absorption cross-section for singlets and triplets. The 
curves obtained with and without the use of symmetry com- 
pletely overlap. 



VI. CONCLUSIONS 



In summary, we presented a scheme that can consid- 
erably reduce the number of calculations necessary to 
obtain the full polarizability tensor by using the sym- 
metries of the system. These can be spatial symmetries 
(like planes of inversion or symmetry axis, e.g.) or they 
can lie in spin space (if the ground-state is spin satu- 
rated). Finally, the scheme is trivial to implement, and 
can be easily extended to different symmetries and dif- 
ferent responses- i.e. higher multipole responses, and 
higher order hyperpolarizabilities -, so we expect that 
its usefulness to surpass the present context. 
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